1. 3D Solid Mechanics#
Open Cascade Techonology geometry kernel
Netgen mesh generator
NGSolve Finite Element libraray
from netgen.occ import *
from netgen.webgui import Draw as DrawGeo
import ngsolve
box = Box((0,0,0), (3,0.6,1))
box.faces.name="outer"
cyl = sum( [Cylinder((0.5+i,0,0.5), Y, 0.25,0.8) for i in range(3)] )
cyl.faces.name="cyl"
geo = box-cyl
DrawGeo(geo);
find edges between box and cylinder, and build chamfers (requires OCC 7.4 or newer):
cylboxedges = geo.faces["outer"].edges * geo.faces["cyl"].edges
cylboxedges.name = "cylbox"
geo = geo.MakeChamfer(cylboxedges, 0.03)
name faces for boundary conditions:
geo.faces.Min(X).name = "fix"
geo.faces.Max(X).name = "force"
DrawGeo(geo);
from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.1)).Curve(3)
Draw (mesh);
1.1. Linear elasticity#
Displacement: \(u : \Omega \rightarrow {\mathbb R}^3\)
Linear strain: $\( \varepsilon(u) := \tfrac{1}{2} ( \nabla u + (\nabla u)^T ) \)$
Stress by Hooke’s law: $\( \sigma = 2 \mu \varepsilon + \lambda \operatorname{tr} \varepsilon I \)$
Equilibrium of forces: $\( \operatorname{div} \sigma = f \)$
Displacement boundary conditions: $\( u = u_D \qquad \text{on} \, \Gamma_D \)$
Traction boundary conditions: $\( \sigma n = g \qquad \text{on} \, \Gamma_N \)$
1.2. Variational formulation:#
Find: \(u \in H^1(\Omega)^3\) such that \(u = u_D\) on \(\Gamma_D\)
holds for all \(v = 0\) on \(\Gamma_D\).
E, nu = 210, 0.2
mu = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))
def Stress(strain):
return 2*mu*strain + lam*Trace(strain)*Id(3)
fes = VectorH1(mesh, order=3, dirichlet="fix")
u,v = fes.TnT()
gfu = GridFunction(fes)
with TaskManager():
a = BilinearForm(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(v))).Compile()*dx)
pre = Preconditioner(a, "bddc")
a.Assemble()
force = CF( (1e-3,0,0) )
f = LinearForm(force*v*ds("force")).Assemble()
from ngsolve.krylovspace import CGSolver
inv = CGSolver(a.mat, pre, printrates='\r', tol=1e-8)
gfu.vec.data = inv * f.vec
CG iteration 1, residual = 0.0001803137821722643
CG iteration 2, residual = 7.486220090884329e-05
CG iteration 3, residual = 8.829232445809833e-05
CG iteration 4, residual = 6.536759069327916e-05
CG iteration 5, residual = 5.845343416703274e-05
CG iteration 6, residual = 4.203085660075135e-05
CG iteration 7, residual = 3.1443770086803845e-05
CG iteration 8, residual = 2.8960104946516772e-05
CG iteration 9, residual = 1.9752033754571705e-05
CG iteration 10, residual = 1.58621922376946e-05
CG iteration 11, residual = 1.0741363250405283e-05
CG iteration 12, residual = 8.121182240847164e-06
CG iteration 13, residual = 6.13958839676602e-06
CG iteration 14, residual = 5.249288611711033e-06
CG iteration 15, residual = 3.791292219643922e-06
CG iteration 16, residual = 2.8087404942947052e-06
CG iteration 17, residual = 2.1122503844054362e-06
CG iteration 18, residual = 1.4856736111195762e-06
CG iteration 19, residual = 1.1232492738871658e-06
CG iteration 20, residual = 7.945075734134604e-07
CG iteration 21, residual = 5.749880111534067e-07
CG iteration 22, residual = 4.1651248011870676e-07
CG iteration 23, residual = 3.149728259459727e-07
CG iteration 24, residual = 2.3300347017451097e-07
CG iteration 25, residual = 1.7753990180064746e-07
CG iteration 26, residual = 1.1969649723751356e-07
CG iteration 27, residual = 9.088435110998614e-08
CG iteration 28, residual = 6.692973749349633e-08
CG iteration 29, residual = 4.8756502489139396e-08
CG iteration 30, residual = 3.5148457873361406e-08
CG iteration 31, residual = 2.5077311983671912e-08
CG iteration 32, residual = 1.7515174615967346e-08
CG iteration 33, residual = 1.3358710552330692e-08
CG iteration 34, residual = 9.017162512655161e-09
CG iteration 35, residual = 6.54138557883829e-09
CG iteration 36, residual = 4.832567168868747e-09
CG iteration 37, residual = 3.6699717291731275e-09
CG iteration 38, residual = 2.4517868480762146e-09
CG iteration 39, residual = 1.925030700827801e-09
CG iteration 40, residual = 1.4219433295295428e-09
CG iteration 41, residual = 9.779828607978307e-10
CG iteration 42, residual = 6.877491752053402e-10
CG iteration 43, residual = 5.264931159192428e-10
CG iteration 44, residual = 3.923016755656033e-10
CG iteration 45, residual = 3.5517945828978116e-10
CG iteration 46, residual = 2.330410842329499e-10
CG iteration 47, residual = 1.5299263810680816e-10
CG iteration 48, residual = 1.2195421670978552e-10
CG iteration 49, residual = 1.0539762592221566e-10
CG iteration 50, residual = 7.757429487426324e-11
CG iteration 51, residual = 5.563405518512371e-11
CG iteration 52, residual = 3.7852906751343235e-11
CG iteration 53, residual = 2.7437899748237933e-11
CG iteration 54, residual = 2.8914277760781623e-11
CG iteration 55, residual = 1.5594882251656798e-11
CG iteration 56, residual = 1.2002867490508141e-11
CG iteration 57, residual = 9.051580388133448e-12
CG iteration 58, residual = 5.900710439121206e-12
CG iteration 59, residual = 4.127934719301097e-12
CG iteration 60, residual = 2.682943853899698e-12
CG iteration 61, residual = 2.05077575758783e-12
CG iteration 62, residual = 1.2837590263734648e-12
CG converged in 62 iterations to residual 1.2837590263734648e-12
with TaskManager():
fesstress = MatrixValued(H1(mesh,order=3), symmetric=True)
gfstress = GridFunction(fesstress)
gfstress.Interpolate (Stress(Sym(Grad(gfu))))
Draw (5e4*gfu, mesh);
Draw (Norm(gfstress), mesh, deformation=1e4*gfu, draw_vol=False, order=3);